################################################################
# R script for replicating conjoint analysis in Legitimacy Book
################################################################


# Load data - US
library(haven)
purpose_us_pooled <- read_dta("~/Desktop/conjointmm/purpose_us_pooled.dta")
View(purpose_us_pooled)

library(rio)
library(gridExtra)
library(ggplot2)
library(cregg)

# conjoint features
tmp <- attr(purpose_us_pooled[["power"]], "label")
purpose_us_pooled[["power"]] <- rio::factorize(purpose_us_pooled[["power"]])
levels(purpose_us_pooled[["power"]]) <- gsub("^.{1}", "???", levels(purpose_us_pooled$power))
attr(purpose_us_pooled[["power"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["purposeR2"]], "label")
purpose_us_pooled[["purposeR2"]] <- rio::factorize(purpose_us_pooled[["purposeR2"]])
attr(purpose_us_pooled[["purposeR2"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["trans"]], "label")
purpose_us_pooled[["trans"]] <- rio::factorize(purpose_us_pooled[["trans"]])
attr(purpose_us_pooled[["trans"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["part"]], "label")
purpose_us_pooled[["part"]] <- rio::factorize(purpose_us_pooled[["part"]])
attr(purpose_us_pooled[["part"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["fair"]], "label")
purpose_us_pooled[["fair"]] <- rio::factorize(purpose_us_pooled[["fair"]])
attr(purpose_us_pooled[["fair"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["problem"]], "label")
purpose_us_pooled[["problem"]] <- rio::factorize(purpose_us_pooled[["problem"]])
attr(purpose_us_pooled[["problem"]], "label") <- tmp
rm(tmp)


# covariates
## internationalist attitudes
tmp <- attr(purpose_us_pooled[["intcoopLOHI"]], "label")
purpose_us_pooled[["Internationalism"]] <- rio::factorize(purpose_us_pooled[["intcoopLOHI"]])
levels(purpose_us_pooled[["Internationalism"]]) <- paste0("Internationalism", levels(purpose_us_pooled[["Internationalism"]]))
attr(purpose_us_pooled[["Internationalism"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["intcoop"]], "label")
purpose_us_pooled[["intcoop"]] <- rio::factorize(purpose_us_pooled[["intcoop"]])
levels(purpose_us_pooled[["intcoop"]]) <- paste0("intcoop: ", levels(purpose_us_pooled[["intcoop"]]))
attr(purpose_us_pooled[["intcoop"]], "label") <- tmp
rm(tmp)

tmp <- attr(purpose_us_pooled[["twolargestparties"]], "label")
purpose_us_pooled[["twolargestparties"]] <- rio::factorize(purpose_us_pooled[["twolargestparties"]])
levels(purpose_us_pooled[["twolargestparties"]]) <- paste0("twolargestparties: ", levels(purpose_us_pooled[["twolargestparties"]]))
attr(purpose_us_pooled[["twolargestparties"]], "label") <- tmp
rm(tmp)

# export
rio::export(purpose_us_pooled, "~/Desktop/purpose_us_pooled.rds")



# calculate quantities for authority and internationalist attitudes

#  Marginal Means
mm <- cj(purpose_us_pooled, confdum ~ power, id = ~ id, weights = ~weight, estimate = "mm", by = ~ intcoop)
mm$statistic <- "MMs"
#  AMCE
amce <- cj(purpose_us_pooled, confdum ~ power, id = ~ id, weights = ~weight, estimate = "amce", by = ~ intcoop)
amce$statistic <- "AMCE"
 

# AMCE results replication - power purpose trans part fair problem
main <- cj(purpose_us_pooled[!is.na(purpose_us_pooled$Internationalism),], confdum ~ power, id = ~ id, weights = ~weight, estimate = "amce_diff", by = ~ Internationalism)
main$statistic <- "Originally Reported Difference"
## fill in baseline difference
main[5L, ] <-  main[4L, ]
main[5L, "level"] <- levels(purpose_us_pooled$power)[1L]
main[5L, c("estimate", "std.error", "z", "p", "lower", "upper")] <- c(0, rep(NA, 5))

# AMCE differences
amce_diffs <- do.call("rbind", lapply(levels(purpose_us_pooled$power), function(x) {
  tmp <- purpose_us_pooled[!is.na(purpose_us_pooled$Internationalism),]
  tmp$power <- relevel(tmp$power, x)
  cj(tmp, confdum ~ power, id = ~ id, weights = ~weight, estimate = "amce_diff", by = ~ Internationalism)
}))
amce_diffs$statistic <- "Potential Differences in AMCEs"

# Difference in Marginal Means
mm_diff <- cj(purpose_us_pooled, confdum ~ power, id = ~ id, weights = ~weight, estimate = "mm_diff", by = ~ Internationalism)
mm_diff$statistic <- "Difference in MMs"

# Merge diffs
diffs <- rbind(main, amce_diffs, mm_diff)
levels(diffs$statistic) <- c("Originally Reported Difference", "Potential Differences in AMCEs", "Difference in MMs")






# calculate quantities for purpose and partisanship

#  Marginal Means
mm <- cj(purpose_us_pooled, confdum ~ purposeR2, id = ~ id, weights = ~weight, estimate = "mm", by = ~ twolargestparties)
mm$statistic <- "MMs"
#  AMCE
amce <- cj(purpose_us_pooled, confdum ~ purposeR2, id = ~ id, weights = ~weight, estimate = "amce", by = ~ twolargestparties)
amce$statistic <- "AMCE"


# AMCE results replication - purposeR2 purpose trans part fair problem
main <- cj(purpose_us_pooled[!is.na(purpose_us_pooled$twolargestparties),], confdum ~ purposeR2, id = ~ id, weights = ~weight, estimate = "amce_diff", by = ~ twolargestparties)
main$statistic <- "Originally Reported Difference"
## fill in baseline difference
main[5L, ] <-  main[4L, ]
main[5L, "level"] <- levels(purpose_us_pooled$purposeR2)[1L]
main[5L, c("estimate", "std.error", "z", "p", "lower", "upper")] <- c(0, rep(NA, 5))

# AMCE differences
amce_diffs <- do.call("rbind", lapply(levels(purpose_us_pooled$purposeR2), function(x) {
  tmp <- purpose_us_pooled[!is.na(purpose_us_pooled$twolargestparties),]
  tmp$purposeR2 <- relevel(tmp$purposeR2, x)
  cj(tmp, confdum ~ purposeR2, id = ~ id, weights = ~weight, estimate = "amce_diff", by = ~ twolargestparties)
}))
amce_diffs$statistic <- "Potential Differences in AMCEs"

# Difference in Marginal Means
mm_diff <- cj(purpose_us_pooled, confdum ~ purposeR2, id = ~ id, weights = ~weight, estimate = "mm_diff", by = ~ twolargestparties)
mm_diff$statistic <- "Difference in MMs"

# Merge diffs
diffs <- rbind(main, amce_diffs, mm_diff)
levels(diffs$statistic) <- c("Originally Reported Difference", "Potential Differences in AMCEs", "Difference in MMs")

